function f=find_T_deriv(x_interest,j,kk,a,b,trunc_low,trunc_high,mu,sigma,fdist,alpha)

N=size(x_interest,1);
f=NaN*ones(N,1);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CASE: (kk=1)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if kk==1

    if j==0
        f(:,1)=0;
    end

    if j>0
        if strcmp(fdist,'Unif')==1
                aux=Tmj_Unif(x_interest,kk,j-1,trunc_low,trunc_high,a,b);
        end
        if strcmp(fdist,'BetaCtr')==1
                aux=Tmj_BetaCtr(x_interest,kk,j-1,alpha,a,b);
        end
        if strcmp(fdist,'TruncNormal')==1
                aux=Tmj_TruncNormal(x_interest,kk,j-1,mu,sigma,trunc_low,trunc_high,a,b);
        end
        if strcmp(fdist,'Normal')==1
                aux=Tmj_Normal(x_interest,kk,j-1,mu,sigma,a,b);
        end

        f(:,1)=j*aux;
    end

end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CASE: (kk=2)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if kk==2
    f(:,1)= -(b+a)/(b-a)* find_T_deriv(x_interest,j,kk-1,a,b,trunc_low,trunc_high,mu,sigma,fdist,alpha) ...
        + 2/(b-a)*find_T_deriv(x_interest,j+1,kk-1,a,b,trunc_low,trunc_high,mu,sigma,fdist,alpha);
end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CASE: (kk>=3)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if kk>=3
    f(:,1)= 4/(b-a)* find_T_deriv(x_interest,j+1,kk-1,a,b,trunc_low,trunc_high,mu,sigma,fdist,alpha) ...
        - (2+4*a/(b-a))* find_T_deriv(x_interest,j,kk-1,a,b,trunc_low,trunc_high,mu,sigma,fdist,alpha) ...
        -  find_T_deriv(x_interest,j,kk-2,a,b,trunc_low,trunc_high,mu,sigma,fdist,alpha);
end

    %
